clear
format short g

% MONOPOLY CASE with entry
% This program looks at the evolution of the computer industry by
% simulating a vintage capital model
%
% version:  (1) Monopolist sells only *1* vintage  
%           (2) each period, products fall (5 to 4, 4 to 3, 3 to 2).
%           (3) Monopolist needs to determine *when* to update product line
%
% story: two kinds of obselecence: (a) outside option (represented by
% falling vintage label and (b) entry of new goods.
%
% Monopolist only chooses when to replace.
% Only one "free" parameter -- fixed cost (which needs to be below the
% total profits number from the competitive case)
%

global inc_flg mean_inc

% Copy output from moments_cmpt.
%int_vec=[b;mc_level; mrkt_size; par_beta1, par_beta2]
par_vec= [8.94861; 0.872581;1.0 ; 0.697482; 1.71831];

%    par_set = [N,nmbr_cmpt,nu_level,a,nu_low,g_rte];
par_set = [10, 15, 10, 1, 0.33, 1.029];

% Income distribution 
   a = par_set(4);
   b = a+par_vec(1);
   par_beta = [par_vec(4); par_vec(5)];
   
% Other parameters
N = par_set(1);
nmbr_cmpt = par_set(2);
nu_level = par_set(3);
nu_low = par_set(5);
g_rte = 1/par_set(6);
  
mc_level = par_vec(2);
mrkt_size = par_vec(3);
% For reporting income stats
    inc_flg = 0;       
    mean_inc = 0;
    
% set of non-frontier number of technologies
    max_rplc=15; 

% discount rate
    par_disc=0.99;

% fixed cost to jumping back to the frontier (make sure it is less than the
% lifetime profits on the competitor's side
% competitors profis: 0.051673
     fx_cost = 0.048;  % 8 periods, matches the data for Apple
% defining vector of vintage utilities
% To create a longer nu_vec, increase N, but then divide the resulting
% nu_vec so that nu_vec(N) = par_set(3).
% defining vector of possible vintage utilities

N = 15; % for price/profit monopolist graph

nu_tmp = zeros(N,1);
nu_tmp(1) = nu_low;
nu_tmp(N) = nu_level;
for i=N-1:-1:2
    nu_tmp(i) = nu_tmp(i+1)*g_rte;
end
nu_vec = nu_tmp;

% Checking to make sure that outside good is least desirable
% If not, then fix
%t_ptr = nu_tmp>par_vec(3);
%nu_vec = [nu_tmp(1);nu_tmp(t_ptr)];
%N = length(nu_vec);

% defining vector of marginal costs
mc_vec = zeros(N,1);
mc_vec(1) = 0;   
mc_vec(2:N) = mc_level;

disp('Full set of utility values and marginal cost');
disp([nu_vec mc_vec]); pause

% Cycling over gaps of 1 through max_rplc-1
rec_outcomes = zeros(max_rplc+1,3);


for k = 1:max_rplc
    
    disp(['Looking at vintage :' num2str(k)]);
    
    nu_t = [nu_vec(1); nu_vec(N+1-k,1)];
    mc_t = [mc_vec(1); mc_vec(N+1-k,1)];

     disp('Utility and marginal cost vectors of potential goods');
     disp([nu_t mc_t]); %pause
    
    % setting initial prices (for one technology, use marginal cost)
    initial_prc = mc_t(2);
    
    % searching for optimal prices
    % a lot of this code is not necessary for the 1 vintage case, but has
    % been taken from code looking at the multiple vintage case
    flg=0;
    iter=1;
    while flg<1
 %       disp(['Scenario & iteration:    ' num2str([ k iter])]);
        prc_equil = fminsearch(@(x) mply_prft(x,nu_t,mc_t,a,b,par_beta), initial_prc);
        
        pp = [0;prc_equil];
        Qd = demand(pp,nu_t,a,b,par_beta);
        Qd_short = Qd(2:length(Qd));
        indx1 = find(Qd_short>0);
        indx2 = find(Qd_short<=0);
        i_pp = [0;initial_prc];
        dd1 = sum(abs(initial_prc(indx1)-prc_equil(indx1)));
        dd2 = sum((pp-mc_t).*Qd);
%        disp(['sum of differences in price and sum of profits:  ' num2str([dd1 dd2])]);
%        disp('Initial and equilibrium price vector');
%        disp([initial_prc prc_equil]);
        if(dd1>0)  % repeat search - make sure zero demand product prices are re-initialized
            ll = length(prc_equil);
            % if no good is selling, reset prices to marginal cost
            if(length(indx1)<1)
                initial_prc = mc_t(2:ll+1);
            else
                tmp_prc = zeros(ll,1);
                tmp_prc(indx1) = prc_equil(indx1);
                % if necessary, filling in highest utility price
                if(tmp_prc(ll)<=0)
                    tmp_prc(ll) = prc_equil(max(indx1));
                end
                % if lowest utility good does not sell, set price equal to
                % marginal cost
                if(tmp_prc(1)<=0)
                    tmp_prc(1) = mc_t(2);
                end
                % for intermediate goods (use prices of active neighbors)
                for i=2:ll-1
                    if(tmp_prc(i)<=0)
                        i_tmp = find(tmp_prc>0);
                        tmp = find(i_tmp<i); low = max(i_tmp(tmp));
                        tmp = find(i_tmp>i);  high = min(i_tmp(tmp));
                        tmp_prc(i) = (tmp_prc(low)+tmp_prc(high))/2;
                    end
                end                
                initial_prc = tmp_prc;
            end
                        
        else
            flg=1;
        end
    end    
    
    % determining static profits
    pp = [0;prc_equil];
    Qd = demand(pp,nu_t,a,b,par_beta);    
    pi_vec = (pp - mc_t).*Qd;
    r = [0; k];
    disp('    Rank    util_val,    price,      mc,    demand,      profit'); 
    disp([r nu_t pp mc_t Qd pi_vec]);

    tmp = sum(pi_vec);
    disp(['Total profits:    ' num2str(tmp)]);
    rec_outcomes(k,1:3)=[tmp prc_equil Qd(2)];
    %pause       
end

disp('vintage Profits prices demand by period');
i = 1:length(rec_outcomes);
disp([i' rec_outcomes]);

% Determining optimal relpacement strategy

V_cycle = zeros(max_rplc+1,2);
V_cycle(:,1) = (1:max_rplc+1)';

flg=0;
k=1;
while flg<1;
%    disp(k);
    j=1;
    for i=0:200
        tmp = rec_outcomes(j,1);
        V_cycle(k,2) = V_cycle(k,2) + (par_disc^i)*tmp;
        if(j>=k)
            j=1;
            V_cycle(k,2) = V_cycle(k,2) - (par_disc^i)*fx_cost;
        else
            j=j+1;
        end
    end
    if(k>1 && V_cycle(k,2)<V_cycle(k-1,2))
        flg=1;
    elseif k>max_rplc
        flg=1; disp('Ran out of vintages');
    else
        k=k+1;
    end
        
end

disp('Replacement Cycle, Total Profits, Prices ');
disp([V_cycle(1:k,:) rec_outcomes(1:k,2)]);

disp('Price and Sales over Product Cycle');
ii = 1:1:(k-1);
disp([ii' rec_outcomes(1:k-1,2:3)]);

    
% Sales moment
sal_cdf = cumsum(rec_outcomes(1:k-1,3)./(sum(rec_outcomes(1:k-1,3))));

data_cdf = [0.112105	0.304142;
            0.304142	0.503962;
            0.503962	0.660601;
            0.660601	0.800873;
            0.800873	0.897171;
            0.897171	0.973031;
            0.973031	1;
            1           1];

            
disp('Data sales CDF with Model sales inbetween'); disp([data_cdf(:,1) sal_cdf data_cdf(:,2)]);

tmp1 = (sal_cdf<=data_cdf(:,1));
tmp2 = (sal_cdf>=data_cdf(:,2));
disp('Number of times outside of the cdf bounds (min, max)'); disp([sum(tmp1) sum(tmp2)]);

disp(' ');
disp(' ');


% Price moment

sal_pnt = 100*sal_cdf;
data_decline = 100*price_declineM(sal_pnt(2:length(sal_pnt)));
p_decline = zeros(1,k-2);
for i=2:k-1
    p_decline(i-1) = 100*(rec_outcomes(i,2)- rec_outcomes(1,2))./rec_outcomes(1,2);
end
disp('Data price decline is'); disp(data_decline');
disp('Model price decline is'); disp(p_decline);

tmp0 = mean(data_decline);
tmp1 = mean(p_decline);
disp(['Mean price declines (unweighted) in data and model:  ' num2str([tmp0 tmp1])]);

disp(' ');
disp(' ');




% Income over the product cycle
sal_rec = rec_outcomes(1:k-1,3);
prc_rec = rec_outcomes(1:k-1,2);

%unlike in the competitive case, only selling one vintage
% so here, looping over the life cycle
mean_inc = zeros(k-1,2);
inc_bnd = zeros(k-1,1);
stp_sze=0.00001;
for i=1:k-1 
    invmshare = 1-sal_rec(i)/mrkt_size;
    mrg_con = icdf('beta',invmshare,par_beta(1),par_beta(2));
    inc_bnd(i) = mrg_con*(b-a)+a;
    mean_inc(i,1) = (inc_bnd(i)+b)/2;
    t1 = mrg_con:stp_sze:1;
    t2 = pdf('beta',t1,par_beta(1),par_beta(2));
    t3 = t1*t2' / sum(t2);
    mean_inc(i,2) = t3*(b-a)+a;    
end

disp('Marginal income by period');
disp([inc_bnd mean_inc]);

